clear;
%% NETWORK-BEHAVIOR COEVOLUTION DATA-GENERATING-PROCESS,
    % and DATA OUTPUT to STATA for CLOGIT MONTE CARLO'S,
    % and DATA OUTPUT to, and ESTIMATION of, SIENA MONTE CARLO'S in R.

    % version: 16 July 2010, 11:35pm, Rob last editor

%% SETUP 
NumTrials=100;    %Number of Monte Carlo Trials
N=50;	%Number of Actors
actors = 1:1:N;	%vector of actor numbers, counting 1 to N by N
T=20;	%Number of Over-Time Observations (beyond the 1st) on the Network (and Behavior)
rate_net = 5;	%network rate of change within intervals ("micro-steps")
rate_beh = 5;	%behavior rate of change within intervals ("micro-steps")
behmat_obs=zeros(N,T+1);	%preallocate Nx(T+1) matrix observed behaviors
netmat_obs=zeros(N,N,T+1);  %preallocate NxNx(T+1) matrix (i.e., three-dimensional array) of observed networks. (Note: SIENA will not accept a sociomatrix with all zeros).
for z = 1:N
	if z <= (N+1)/2
        behmat_obs(z,1) = 1;	%set first .5N actors' behaviors in first-period Nx1 vector of observed behaviors to 1
        if z < N/2
            netmat_obs(z,z+1)=1;	%each of these first .5N actors are connected to next actor, except last... who connected to 1st
        else
            netmat_obs(z,1)=1;  %...except last, who connected to 1st
        end
    else
        behmat_obs(z,1) = 0;	%set second .5N actors' behaviors in first-period Nx1 vector of observed behaviors to 0
        if z < N
            netmat_obs(z,z+1)=1;	%each of these second .5N actors are connected to next actor, except last...
        else
            netmat_obs(z,1)=1;  % ...except last, who connected to 1st
        end
	end
end
m = zeros(T+1,1);   %initializing m, which will record the number of the micro-steps at which the observed periods arise
m(1,1)=1;
SL = zeros(N,T);    %initializing SL, the spatial lag, to an NxT matrix.
results_logit_beh = zeros(NumTrials,6);   %initialize simple-estimator behavior-equation output to (NumTrials x 6); note: 6=(3 coefficients in behavior eqtn * 2 statistics, coeff & se est)
results_logit_tie = zeros(NumTrials,6);   %initialize simple-estimator network-tie-equation output to (NumTrials x 6); note: 6=(3 coefficients in tie eqtn * 2 statistics, coeff & se est)

%%  Big Loop for Montes 
for zz = 1:NumTrials   %zz counts the trials
    clear ch_combT; %clearing ch_combT from previous trial;
    ch_combT = [0 0 0 0 0]; %initializing ch_combT; (eventually) a matrix storing, by column:
                            %col1: micro-steps at which changes made,
                            %col2: whether change to network (0) or behavior (1),
                            %col3: actor number selected to make change,
                            %col4: network-tie number changed by actor (+ = switched on; - = switched off),
                            %col5: behavior-change made by actor (0=no change; +1=changed to on; -1= changed to off).
%% Determine who changes and when; See Snijders (2001), pp 367-368.
for t=1:T
    ch_t_net(1) = 0; %initialize first-period (micro-step) probability of a network-change event to zero
    ch_t_beh(1) = 0; %initialize first-period (micro-step) probability of a behavior-change event to zero
    i=1;
    while (ch_t_net(i) < 1)
        ch_t_net(i+1) = ch_t_net(i) + (-log(1-rand))/(N*rate_net); %create vector of cumulative probabilities network-change event (micro-steps), following an exponential distribution
        i = i+1;
    end
    ch_t_net_comb = zeros(2,i-2); %preallocate 2x(i-2) placeholder for matrix of event-probabilities (last two probabilities(i) generated by preceding loop exceed 1)
    ch_t_net_comb(1,:) = ch_t_net(2:i-1); %first row of that event-probability matrix are the network-change event-probabilities just generated (retaining zero for first-step).
    i=1;
    while (ch_t_beh(i) < 1)
        ch_t_beh(i+1) = ch_t_beh(i) + (-log(1-rand))/(N*rate_beh); %create vector of cumulative probabilities behavior-change event (micro-steps), following an exponential distribution
        i = i+1;
    end
    ch_t_beh_comb = ones(2,i-2); %preallocate 2x(i-2) placeholder for matrix of event-probabilities (last two probabilities(i) generated by preceding loop exceed 1)
    ch_t_beh_comb(1,:) = ch_t_beh(2:i-1); %first row of that event-probability matrix are the behavior-change event-probabilities just generated (retaining zero for first-step).
    ch_comb = [ch_t_net_comb, ch_t_beh_comb]'; %concatenate the network-change and the behavior-change event-probability vectors, transposing to (number-of-steps)x2 matrix
    ch_comb(:,3) = randsample(actors, length(ch_comb), true); % add third column to this matrix that is a uniform random draw, with replacement, from the N actors (see Snijders 2001 top p. 368)
    ch_comb(:,4) = zeros(length(ch_comb),1);
    ch_comb(:,5) = zeros(length(ch_comb),1);
    ch_comb = sortrows(ch_comb,1);	%ch_comb sorted by its first column, which is micro-step, measured in percent of way from last observation to next, of events; yields the following record of unobserved changes:
                                    %[micro-steps at which changes made, whether change to network (0) or behavior (1), and actor number selected to make change]
    m(t+1,1) = m(t,1)+length(ch_comb);  %m keeps record of number of micro-steps at which next observation occurs
    ch_combT = [ch_combT; ch_comb]; %vertical concatenation of ch_comb from the intraobservation periods
end
ch_combT = ch_combT(2:length(ch_combT),:);  %drop first row of ch_combT (which all zeros)
behmat_total=zeros(N,length(ch_combT)+1);	%preallocate matrix of behaviors, dimensions: Nx(number of microsteps)
netmat_total=zeros(N,N,length(ch_combT)+1);	%preallocate array of networks, dimensions: NxNx(number of microsteps)
behmat_total(1:N,1) = behmat_obs(1:N,1);	%first column of behaviors, at microstep 1, are the first-observed behaviors
netmat_total(1:N,1:N,1) = netmat_obs(1:N,1:N,1);	%first network matrix, at microstep 1, is the first-observed network

%% Generate Networks and Behavior
for t=1:T

%Determine the changes using the Objective Functions and Random Utility Model 
%Calculate Change Probabilities for Network Ties and Behavior using Objective Functions (Currently Implementing Covariate-Related Similarity and Average-Similarity, see RSiena manual pp. 56 and 60)
	for i=m(t):m(t+1)-1     %Loop for Unobserved Data, i.e. over all microsteps
    %CODE FOR ACTORS CHANGING THEIR NETWORK TIES:
        if ch_combT(i,2) == 0	%if ch_comb 2nd col for this microstep is 0, then event is a network-change possibility. (Actually, change will be made as coded, following Snijders (2001).)
            beh_range = range(behmat_total(1:N,i));	%finds range, max minus min, of the N actors' behaviors at step i
            diff_mat=zeros(N,N);    %preallocate NxN matrix zeros for differences in actors' behaviors at step i
            for j=1:N
                diff_mat(1:N,j) = behmat_total(1:N,i) - behmat_total(j,i);  %calculate differences in actors' behaviors at step i
            end
            if beh_range == 0   %if all actors same behavior in step i
                sim_mat = ones(N,N);    %then similarity matrix all ones
                av_sim = 1; %and average similarity is one
            else
                sim_mat =(beh_range - abs(diff_mat))/beh_range; %else elements of similarity matrix are 0..1 with 0 if maximally, for step 1, different and 1 if same
                av_sim = (sum(sum(sim_mat))-N*beh_range)/(N*(N-1)); %average similarity score across all actors.
                %UNCLEAR: average for actor i or average over all actors? (see manual, p.56)
            end
            stat_net_n=zeros(1,N);    %preallocate 1xN vector zeros for actor j's network-statistics (for covariate-similarity selection) at step i
            for j=1:N 
                if j==ch_combT(i,3);    %for the actor chosen to make network-change at step i...  
                    stat_net_n(j) = 0;  %...leave network-stat zero.
                else    % Implement Covariate-Related Similarity (see RSiena manual pp. 56 and 60)
                    netmat_total_temp = netmat_total(ch_combT(i,3),:,i);    %For the jth actor's set of covariate-similarity scores...
                    netmat_total_temp(1,j,i) = 1 - netmat_total(ch_combT(i,3),j,i); %...take the dissimilarity in their behaviors and...
                    stat_net_n(j) = netmat_total_temp(1,:,i)*(sim_mat(ch_combT(i,3),:)'-av_sim); %...calculate the relative covariate-similarity scores of j with all others.
                end       
            end
            probs_ru=zeros(1,N);    %preallocate 1xN vector zeros for actor j's probabilities of changing each network connection at step i
            for j=1:N
                if j~=ch_combT(i,3);    %For each actor not j...
                    probs_ru(j) = exp(stat_net_n(j))/(sum(exp(stat_net_n))-exp(stat_net_n(ch_combT(i,3))));     %Calculate the (Multinomial Logit) Probabilities of j Changes each Network Tie.
                                                                                                                %NOTE: Coefficient of 1 on the network-statistic hardwired here.
                else
                    probs_ru(j)= 0;     %Self-ties are not accomodated in this model.
                end
            end 
            net_options_n = 1:1:N;  %Initialize another counter vector, 1 to the number of actors.
            p = probs_ru/sum(probs_ru);     %Fix for 2009b randsample.m bug, added by Rob 14 July 2010
            e = min([0 cumsum(p)],1);     %Fix for 2009b randsample.m bug, added by Rob 14 July 2010
            e(end) = 1;     %Fix for 2009b randsample.m bug, added by Rob 14 July 2010
            p = diff(e);     %Fix for 2009b randsample.m bug, added by Rob 14 July 2010
            ch_net = randsample(net_options_n, 1, true, p);         %Choice of Network Tie Change -- a random draw, with probabilities given by MNL model above, of the tie j will change.
            netmat_total(1:N,1:N,i+1) = netmat_total(1:N,1:N,i);	%Network for microstep i+1 is Network from microstep i, except...
            netmat_total(ch_combT(i,3),ch_net,i+1) = 1 - netmat_total(ch_combT(i,3),ch_net,i);	%...for the change j just made in its network ties.
            behmat_total(1:N,i+1) = behmat_total(1:N, i);           %Behavior vector for microstep i+1 remains unchanged at behavior vector from i.
            ch_combT(i,4) = (netmat_total(ch_combT(i,3),ch_net,i+1)-netmat_total(ch_combT(i,3),ch_net,i))*ch_net;   %Stores + or - or 0 changed tie just made in column 4 of ch_combT.
        else	%if ch_comb 2nd col for this microstep was 0, then event is a behavior-change possibility
    %CODE FOR CHANGE IN BEHAVIOR
            beh_range = range(behmat_total(1:N,i));	%finds range, max minus min, of the N actors' behaviors at step i                                       
            for j=1:N
                diff_mat(1:N,j) = behmat_total(1:N,i) - behmat_total(j,i);  %calculate differences in actors' behaviors at step i
            end
            if beh_range == 0   %if all actors same behavior in step i
                sim_mat = ones(N,N);    %then similarity matrix all ones
                av_sim = 1; %and average similarity is one
            else
                sim_mat =(beh_range - abs(diff_mat))/beh_range; %else elements of similarity matrix are 0..1 with 0 if maximally, for step 1, different and 1 if same
                av_sim = (sum(sum(sim_mat))-N*beh_range)/(N*(N-1)); %average similarity score across all actors.
                %UNCLEAR: average for actor i or average over all actors? (see manual, p.56)
            end
        %Utility from No Change in Behavior
            if sum(netmat_total(ch_combT(i,3),:,i))==0	%If actor chosen to consider change behavior is unconnected to anyone...
                stat_beh_nc_n=0;    %...then the utility from no-change is 0,
            else
                stat_beh_nc_n = (netmat_total(ch_combT(i,3),:,i)*(sim_mat(ch_combT(i,3),:)'-av_sim))/(sum(netmat_total(ch_combT(i,3),:,i)));    %...else utility of no-change is covariate-similarity at microstep i.
            end
        %Utility from Change in Behavior
            behmat_total_temp = behmat_total(1:N,i);    %Create a hypothetical vector behaviors equal to step i's behaviors except that...
            behmat_total_temp(ch_combT(i,3),1) = 1 - behmat_total(ch_combT(i,3),i); %...the chosen actor's behavior is switched.
            beh_range_temp = range(behmat_total_temp(1:N,1));   %Calculate the behavior range in this hypothetical behavior vector.
            diff_mat_temp=zeros(N,N);   %Preallocate an NxN matrix for the differences in the network from this behavior change.
            for j=1:N
                diff_mat_temp(1:N,j) = behmat_total_temp(1:N,1) - behmat_total_temp(j,1);   %Calculate that matrix of behavioral differences.
            end
            if beh_range_temp == 0  %If this hypothetical behavior-change aligns all N behaviors...
                sim_mat_temp = ones(N,N);   %...then the similarity matrix is ones...
                av_sim_temp = 1;    %...and the average similarity is one.
            else
                sim_mat_temp =(beh_range_temp - abs(diff_mat_temp))/beh_range_temp; %Otherwise, calculate the behavioral-similarity matrix under this hypothetical change.
                av_sim_temp = (sum(sum(sim_mat_temp))-N*beh_range_temp)/(N*(N-1));  %Calculate the average similarity score under this hypothetical change.
                                                                                    %UNCLEAR: average for actor i or average over all actors? (see manual, p.56)
            end
            if sum(netmat_total(ch_combT(i,3),:,i))==0	%If actor chosen to consider change behavior is unconnected to anyone...
                stat_beh_ch_n=0;	%...then the utility from no-change is 0,
            else
                stat_beh_ch_n = (netmat_total(ch_combT(i,3),:,i)*(sim_mat_temp(ch_combT(i,3),:)'-av_sim_temp))/(sum(netmat_total(ch_combT(i,3),:,i)));
                %...else utility of change is covariate-similarity under hypothetical change.
            end
            stat_beh_n = [stat_beh_nc_n, stat_beh_ch_n];    %Store the utilities to compare, of no-change & change, in a vector.
            probs_beh = exp(stat_beh_n)/(sum(exp(stat_beh_n))); %(Multinomial, but just 2 options, so) Logit Random-Utility Model of Probability of Behavior Change.
                                                                %NOTE: Coefficient of 1 on the network-statistic hardwired here.
            beh_options_n = 1:1:2;  %Create a vector of the possible behaviors, 1 or 2.
            ch_beh = randsample(beh_options_n, 1, true, probs_beh); %Randomly choose the behavior based on (Logit) Random-Utility Model.
            if ch_beh == 1
                behmat_total(1:N,i+1) = behmat_total(1:N,i);	%If j not change behavior, then behavior vector i+1 = vector i.
            else
                behmat_total(1:N,i+1) = behmat_total(1:N,i);	%If j change behaviors, then behavior vector i+1 = vector i, except...
                behmat_total(ch_combT(i,3),i+1) = 1 - behmat_total(ch_combT(i,3),i);    %...for row j where the change made.
            end
            netmat_total(1:N,1:N,i+1) = netmat_total(1:N,1:N,i);    %The network remains the same from i to i+1 because this event was a behavior-change possibility.
            ch_combT(i,5) = behmat_total(ch_combT(i,3),i+1) - behmat_total(ch_combT(i,3),i);   %Stores +1 or -1 or 0 for changed behavior just made in column 5 of ch_combT. 
        end
	end %END of i LOOP (i.e, the loop for Unobserved DATA, the microsteps)
    %Once reach end of the m(t)-to-m(t+1) microsteps...
    behmat_obs(1:N,t+1) = behmat_total(1:N,i+1);    %...record behavior vector at that microstep as the OBSERVED BEHAVIORS for t+1.
    netmat_obs(1:N,1:N,t+1) = netmat_total(1:N,1:N,i+1);    %...and record network matrix at that microstep as the OBSERVED NETWORK for t+1.
    SL(1:N,t) = normw(netmat_obs(:,:,t))*behmat_obs(1:N,t); %create spatial-lag of behavior: row-normalized observed-network matrix multiplied by observed-behavior vector.
                                                            %NOTE: SL effectively time-lagged, since for t not t+1
end %END OF BIG-T LOOP (i.e., the loop for Observed DATA)

%% OUTPUT DATA to R for RSiena ESTIMATION
% Data formatting for R (csv)
path_name = 'C:\Users\franzese\Documents\Work\WPDOCS\Spatial Methods\Path Dep\PolMeth 2010\RData\';
for t=1:T+1
    filename = sprintf('%s%s%d%s%d.csv', path_name, 'net_t', t, '_mc', zz);
    csvwrite(filename,netmat_obs(:,:,t));
    filename = sprintf('%s%s%d%s%d.csv', path_name, 'beh_t', t, '_mc', zz);
    csvwrite(filename,behmat_obs(:,t));
end


%% OUTPUT DATA TO, & ESTIMATE SIENA IN, R
%{
% Connect to an R Session
openR

% Load Data into MATLAB and Push into R
putRdata('netmat_obs1',netmat_obs(:,:,1));
putRdata('netmat_obs2',netmat_obs(:,:,2));
putRdata('netmat_obs3',netmat_obs(:,:,3));
putRdata('netmat_obs4',netmat_obs(:,:,4));
putRdata('behmat_obs',behmat_obs)

% Run a series of commands and grab the results
evalR('setwd(''C:/Documents and Settings/Administrator/Desktop/POLMETH_10/MATLAB_R/Results'')');
evalR('library(RSiena)');
evalR('network <- sienaNet(array(c(netmat_obs1,netmat_obs2,netmat_obs3,netmat_obs4),dim=c(30,30,4)))');
evalR('behavior <- sienaNet(behmat_obs, type="behavior")');
evalR('mydata <- sienaDataCreate(network,behavior)');
evalR('myeff <- getEffects(mydata)');
%evalR('fix(myeff)');
evalR('myeff[,9] <- FALSE');
evalR('myeff$include[myeff$functionName%in%"Amount of network change in period 1"] <- TRUE');
evalR('myeff$include[myeff$functionName%in%"Amount of network change in period 2"] <- TRUE');
evalR('myeff$include[myeff$functionName%in%"Amount of network change in period 3"] <- TRUE');
evalR('myeff$include[myeff$functionName%in%"Similarity on behavior"] <- TRUE');
evalR('myeff$include[myeff$functionName%in%"Amount of behavioral change in period 1 on behavior"] <- TRUE');
evalR('myeff$include[myeff$functionName%in%"Amount of behavioral change in period 2 on behavior"] <- TRUE');
evalR('myeff$include[myeff$functionName%in%"Amount of behavioral change in period 3 on behavior"] <- TRUE');
evalR('myeff$include[myeff$type%in%"eval"&myeff$functionName%in%"beh. behavior average similarity"] <- TRUE');
%evalR('fix(myeff)');
evalR('mymodel <- sienaModelCreate(useStdInits = TRUE, projname=''network_evolution'')');
evalR('ans1 <- siena07(mymodel, data=mydata, effects=myeff)');
theta = evalR('ans1$theta');
covtheta = evalR('ans1$covtheta');
results(zz,:) = [theta,sqrt(diag(covtheta))'];


% Close the connection
closeR
%}

%% Reorganize DATA and estimate simple logits
%BEHAVIOR MODEL
beh_logit = reshape(behmat_obs, N*(T+1), 1);    %Reshape Nx(T+1) observed-behavior matrix to N(T+1)x1 vector
beh_logit_templag = beh_logit(1:(N*(T+1)-N), 1);    %Create 1-observation-period lag of observed-behavior vector
beh_logit = beh_logit(N+1:N*(T+1), 1);	%Drop the first observation-period (for which do not have lag & SIENA model cannot be applied)
beh_logit_TLSL = reshape(SL, N*T, 1);   %Reshape (time-lagged) spatial-lag to NTx1.

X1 = [ones(N*T,1),beh_logit_templag, beh_logit_TLSL];   %Regressors of time-lagged spatial-lag logit behavior-model

result1 = logit(beh_logit,X1);  %Estimate time-lagged spatial-lag logit behavior-model.
tstat1 = result1.tstat; %Save t-stats
beta1 = result1.beta;   %Save coefficient estimates
se1 = beta1./tstat1;    %Save standard-error estimates
results_logit_beh(zz,:) = [beta1', se1'];   %Enter results for this trial (in order: b_intercept, b_timelag, b_spatiallag, se_intcpt, se_tlag, se_splag)
% Save Logit-Behavior Results (csv)
path_name = 'C:\Users\franzese\Documents\Work\WPDOCS\Spatial Methods\Path Dep\PolMeth 2010\RData\';
filename = sprintf('%s%s%d%s%d%s%d%s%d%s%d.csv', path_name, 'logit_beh_n', N, '_t', T, '_Rb', rate_beh, '_Rn', rate_net, '_mc', zz);
    csvwrite(filename,results_logit_beh(zz,:));

%NETWORK MODEL

%Reshape NxNx(T+1) array of observed networks to (N(N-1))x(T+1) matrix
dyad=zeros(N*(N-1),T+1);	%initialize matrix to hold the (N(N-1))x(T+1) matrix
for t=1:T+1
    k=1;
for i = 1:N
    for j = 1:N
        if i~=j
        dyad(k,t)=netmat_obs(i,j,t);
        k=k+1;
        else
        end
    end
end
end

homophily=zeros(N*(N-1),T+1);	%create (N(N-1))x(T+1) 'homophily matrix' with element ij=1 if behavior_i=behavior_j, =0 else
for t=1:T+1
    k=1;
for i = 1:N
    for j = 1:N
        if i~=j
        if behmat_obs(i,t)==behmat_obs(j,t); 
                homophily(k,t) = 1;
                k=k+1;
        else
                homophily(k,t) = 0; 
                k=k+1;
        end
        else
        end
    end
end
end

%Reshape these N(N-1)x(T+1) matrices to N(N-1)(T+1)x1 vectors
%Create time-lag of the ties.
%Save time-lagged 'homophily vector'
tie_logit = reshape(dyad(1:(N*(N-1)),2:T+1),N*(N-1)*T,1);
tie_logit_templag = reshape(dyad(1:(N*(N-1)),1:T),N*(N-1)*T,1);
homophily_logit = reshape(homophily(1:(N*(N-1)),1:T),N*(N-1)*T,1);

X2 = [ones(N*(N-1)*T,1), tie_logit_templag, homophily_logit];   %Regressors of time-lagged-homophily network-tie model
result2 = logit(tie_logit, X2);
tstat2 = result2.tstat;
beta2 = result2.beta;
se2 = beta2./tstat2;
results_logit_tie(zz,:) = [beta2', se2'];
% Save Logit-Tie Results (csv)
path_name = 'C:\Users\franzese\Documents\Work\WPDOCS\Spatial Methods\Path Dep\PolMeth 2010\RData\';
filename = sprintf('%s%s%d%s%d%s%d%s%d%s%d.csv', path_name, 'logit_tie_n', N, '_t', T, '_Rb', rate_beh, '_Rn', rate_net, '_mc', zz);
    csvwrite(filename,results_logit_tie(zz,:));
%% END of MONTE CARLO LOOP
trial = zz
end